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ABSTRACT 

We have constructed a new code to produce synthetic spectra of stellar populations 
that includes massive binaries. We have tested this code against the broadband colours 
of unresolved young massive stellar clusters in nearby galaxies, the equivalent widths 
of the Red and Blue Wolf-Rayet bumps in star- forming SDSS galaxies and the UV and 
optical spectra of the star forming regions Tol-A and B in NGC5398. In each case we 
find a good agreement between our models and observations. We find that in general 
binary populations are bluer and have fewer red supergiants, and thus significantly 
less flux in the I-band and at longer wavelengths, than single star populations. Also we 
find that Wolf-Rayet stars occur over a wider range of ages up to 10^ years in a stellar 
population including binaries, increasing the UV flux and Wolf-Rayet spectral features 
at later times. In addition we find that nebula emission contributes significantly to 
these observed properties and must be considered when comparing stellar models 
with observations of unresolved stellar populations. We conclude that incorporation of 
massive stellar binaries can improve the agreement between observations and synthetic 
spectral synthesis codes, particularly for systems with young stellar populations. 

Key words: galaxies: starburst - galaxies: star clusters - galaxies: stellar content - 
binaries: general - stars: evolution - stars: Wolf-Rayet 



1 INTRODUCTION 

While bright individual stars can be resolved by observa- 
tions of local galaxies, for more compact or distant systems 
only the ensemble properties of an unresolved population 
can be measured. Through necessity, the analysis applied 
in these cases differs. For resolved populations it is possible 
to study each star in detail, and to compare number ratios 
of different stellar types. By contrast, in unresolved popu- 
lations the relative contributions of different stellar types 
are estimated from the overall spectral energy distribution 
and from specific emission/absorption lines by using syn- 
thetic colours and spectra calculated from spectral popula- 
tion synthesis (SPS) codes. However SPS codes suffer from 
many uncertainties due to limited modelling of contributions 
from numerically small, but nonetheless important stellar 
sub- types, and due to different assumptions and techniques 
for balancing t he contribution of different c omponents in 
the population l|Conrov. Gunn fc Whitell2008l ). Binary evo- 
lution, and in particular the effects of binary evolution on 
massive, short-lived but ultraviolet luminous stars such as 
the Wolf-Rayet population, has been neglected in existing 
SPS codes. 
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An understanding of the signatures of massive stars in 
unresolved young populations is essential for correctly in- 
terpreting the star formation hist ory and environment of 
starburst regions both lo cally (sec Sc haerer fc Vaccalll998l : 
iBrinchmann. Kunth fc Durret 200 8) and at high redshifts 
where starbursts are more prevalent, and data - both photo- 
metric and spectroscopic - sparser. At the highest redshifts, 
the predicted s ignatures of massive Population III stars (e.g. 
ISchaereil|20 02!) must be disentangled from those of massive 
stars at higher metallicities, and so the latter must first be 
well understood. In this paper we address the inclusion of 
binary evolution paths in spectral synthesis codes, and ex- 
plore their effects on fitting of unresolved stellar spectra of 
star-forming galaxies. 

lEldridge, Izzard fc ToutI (|2008l ) demonstrated that the 
predicted properties of stellar populations including bi- 
nary models matched those of resolved stellar pop- 
ulations better than single-star m odels. Independently 
IBrinchmann. Kunth fc Purred l|2008l ) showed that the same 
binary models provide a good fit to the stellar populations of 
unresolved Wolf-Rayet galaxies in Sloan Digital Sky Survey 
(SDSS) data. Their analysis derived the relative numbers of 
O stars and Wolf-Rayet stars from the equivalent widths of 
certain spectral features in the optical regime, and compared 
this ratio with those modelled. A more direct comparison is 
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possible: between the observed spectra and synthetic spectra 
derived from a model binary population. 

Codes to create such synt hetic spectra ex ist, the most 
widely-used being starburst99 (jLeitherer et al.lil999, ). How- 
ever few such codes take into consideration the full effect 
of binary ev o lution . IVan Bever fc VanbeverenI (Il998l ) and 
iBelkus et al] (|2003l ) and references therein included binary 
evolution in a rapid population synthesis code. Their bi- 
nary model was based o n 1000 detailed stellar models, while 
iDionne fc Ro'totI l|2006D adapted the same binary models to 
be included in starhurst99. However since the single star and 
binary models came from different stellar evolution codes 
they modified the binary models to match with the sin- 
gle star evolut i onary tracks in starburst99. We also note 
that iHan et all l|2007l ) incorporated binary evolution of low- 
mass stars to explain the observed UV flux in some ellipti- 
cal galaxies. However they have not extended this to higher 
mass binaries. 

Our self-consistent approach uses full sets of detailed 
single stars and binary models from lEldridge. Izzard fc Tout! 
(|2008 l). Both were created with one code, the Cambridge 
STARS code, with identical input physics. Also the evo- 
lution of our binary population is based on 3000 detailed 
models for each metallicity. 

In addition our spectral synthesis code presented here 
incorporates theoretical rather than empirical stellar spec- 
tral libraries wherever possible in order to provide as com- 
plete a theoretical model as possible to compare to obser- 
vational data. Therefore the only empirical results that go 
into our code are the equivalent widths of HeH lines for Of 
stars and a subset of WNL stars, the strength of convec- 
tive overshooting, the mixing length for convection and the 
mass-loss rates for red supergiants and Wolf-Rayet stars. 

Finally, given this theoretical bias, we account for neb- 
ular emission from the gas and dust surrounding our stel- 
lar popula tion through use of the photoionization program 
Cloudy (|Ferland et al.1 Il998l ). In combination with the 
purely synthetic stellar spectra, this leads to an accurate 
model of the total emission from stars, gas and dust in a 
model galaxy or star cluster. 

The structure of this paper is as follows: In section [5] 
we detail the construction of our synthetic stellar spectra, 
and how we process this output with Cloudy. In section 
Owe compare the broad-band magnitudes and colours pre- 
dicted from our models to observed young massive clusters 
in nearby galaxies (section I3.1|l , to equivalent line widths 
for Wolf-Rayet features in unresolved SDSS galaxies (sec- 
tion |3]2]), and to the observed UV and optical spectra of the 
Tol A fc B regions discussed by ISidoli. Smith fc Crowthen 
l|2006l , section [3^ . Finally we discuss the implications and 
interpretation of our results and list our conclusions. 



2 SYNTHETIC SPECTRA OF STELLAR 
POPULATIONS 

Creating a synthetic spectra of a stellar population requires 
the combination of stellar evolution models, model stellar 
atmospheres and a nebular emission model. In this section 
we detail the source of these sets of models and how they 
are combined to produce a composite population spectrum. 



2.1 Stellar evolution models 

We us e stellar models from the Cambridge S TARS 
code (|Eggletonl Il97ll : lEldridge. Izzard fc Tou3 |2008| . 
and r eferences therein), spec ifically those calculated in 
lEldrid go. Izzard & Tout (|2008h . Their key feature is that 
there is not only a set of detailed single star models, but 
also an extensive set of detailed binary star models which 
are key to producing a realistic synthetic stellar population. 
We consider stellar models at five different metallicities: 
Z^O.OOl, 0.004, 0.008, 0.020 and 0.040 (where a metallicity 
of Z=0.020 is conventionally considered solar), with hy- 
drogen mass fraction, X = 0.75 — 2.5Z, and helium mass 
fraction, Y = 0.25 -I- 1.5Z. We use the method described in 
lEldridge. Izzard fc ToutI { 20081 ) to model the primary and 
secondary stars in a binary and to account for mergers. 
The only difference in our current procedure is that at 
each timestep we now select a model atmosphere that best 
represents the model and combine these to form a composite 
spectrum for the population. Given that stellar evolution 
is non-linear and binary evolution is even less predictable, 
we do not interpolate between models with different masses 
and initial binary parameters, but rather weight each 
stellar model by a Salpeter initial mass function (IMF). We 
note t hat the results presented in lEldridge. Izzard fc Tout! 
(|2008t ) are for continuous star formation. Here we consider 
the evolution of a single instantaneous burst with age, 
and do not consider multiple bursts of star formation. 
Features in the rest-UV are sensitive primarily to young 
stars so are relatively unaffected by older underlying stellar 
populations, while such an older stellar population will 
tend to boost the optical continuum, and hence reduce the 
strength of line emission features in this region. 



2.1.1 Description of the binary models and population 
synthesis 

While a full description of our d etails models can be found in 
lEldridge. Izzard fc ToutI (2008) we provide a brief overview. 
We have modified our stellar evolution code to model bi- 
nary evolution. The details of our binary interaction algo- 
rithm are relatively simple compared to the sc heme out- 
lined in, for example. iHurlev. Tout fc Fold l|2002l ). We make 
a number of assumptions in producing our code to keep it 
relatively simple. Our aim was to investigate the effect of 
enhanced mass loss due to binary interactions on stellar life- 
times and populations; therefore we concentrated on these 
aspects rather than incorporating additional physical pro- 
cesses, each of which would add more free parameters to our 
models and potentially associated uncertainties on those pa- 
rameters or the mechanisms concerned. For example we do 
not include wind accretion, gravitational breaking or mag- 
netic breaking. The processes are unlikely to be important 
in the evolution of massive binaries due to the short evolu- 
tionary timescales of massive stars. Processes that would 
have a more significant results on our binary population 
are that we employ a simple tidal model that assumes stel- 
lar rotation only becomes synchronous with the orbit dur- 
ing mass-transfer events and t hat all the orbits are circu- 
lar t hroughout their evolution. iHurlev. Tout fc Poi3 (|2002l ) 
and IStancliffe fc Eldridgd l|2009l ) find that including accu- 
rate treatments of these do not provide new evolutionary 
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paths, but do alter the initial separation at which different 
evolutionary paths occur. For example inclusion of an accu- 
rate tide model increases the maximum initial separation for 
a mass-transfer event to occur in a binary. This suggests we 
could be underestimating the number of binary interactions 
in our stellar models. 

We also make assumptions in calculating our synthetic 
population to avoid calculating a large number of models. 
For example, we do not model accretion onto the secondary 
in the detailed code. We take the greater of the final mass 
of the secondary at the end of the primary code or its initial 
mass when we consider the evolution of the secondary. This 
avoids calculating 10 times more secondary models than pri- 
mary models but possibly missed important effects of accre- 
tion on the evolution. 

We always treat the primary as the initially more mas- 
sive star and we only evolve one star at a time with our 
detailed code. Our models have initial separations that take 
values between log(o/-R0) = 1 to 4 in steps of 0.25 dex. The 
mass ratio takes values of q — 0.1, 0.3, 0.5, 0.7 and 0.9. The 
primary initial masses range from 5 to 120 M0, so the small- 
est mass star in our population is 0.5 Mq. However we only 
include a star's contribution to our synthetic spectrum if its 
mass is above 5 Mq at any point of its evolution. Because 
of this constraint, we do not attempt to model populations 
older than approximately 40 Myrs old. 

When we evolve the primary in detail, it has a shorter 
evolutionary time-scale than the secondary which remains 
on the main sequence until after the primary completes 
its evolution and so we can determine the state of the 
secondary using the singl e stellar evolution equations of 
iHurlev. Pols fc ToutI (|2000l ). When we evolve the secondary 
in detail, we assume that its companion is the compact rem- 
nant of the primary (a white dwarf, neutron star or black 
hole) and treat this as a point mass. If we find that the 
binary systems experience a merger we use a very simple 
merger scheme in which, when the stars come into contact, 
all the mass of the secondary is accreted onto the primary. 
Then subsequent evolution occurs as a single star. 

Our population synthesis calculations a re built upon 
those described in lEldridge. Izzard fc ToutI l|2008l ') with 
some improvements. The important point of evolution is 
what happens after the first supernova (SN) in a binary. If 
a star has a carbon/oxygen core mass greater than 1.38 M© 
and the final mass of the star is greater than 2Mq we 
assume it explodes in a SN. If a neutron star is formed, 
we determine the fate o f the b inary by using the work of 
iBrandt fc PodsiadlowskH (|l995l ) with the latest determina- 
tion for the k ick ve locity distribution from observations by 
iHobbs et all (|2005l ). If the remnant is a black hole, we as- 
sume that it receives a similar kick. Because the masses of 
black holes are greater than neutron stars we use the kick 
distribution of Hobbs et al. ( 2005) to fix the momentum dis- 
tribution. We calculate the resulting black hole kick velocity 
from VBH = ■yNs(l-4/MBH). 

For each primary model there are many possible out- 
comes after the first SN. The binary might be unbound or 
remain bound. In the latter case there are a range of differ- 
ent orbital separations possible, dependant on the strength 
and direction of the SN kick. To determine how important 
the different possible outcomes are we generate a random 
SN kick velocity and direction, calculate the effect on the 



binary, then repeat many times to estimate the relative im- 
portance of the outcomes. This process leaves us with the 
weights to apply to our secondary models when we include 
them in our synthetic population and spectra. 

2.1.2 The effect of rotation 

Modelling binary systems provides in general similar effects 
to including rotation in single star models. Both are very 
complex to implement and currently there are few codes that 
include both. To our knowl edge only the c ode described by 
ICantieUo et~aLl (|2007| ) and Ide Mink etHI 12009,) (and ref- 
erences there-in) does so, and is extremely computationally 
and human input intensive. Trying to separate out the effects 
of binaries and single-star rotation on stellar populations 
is difficult as both have similar effects, that is to produce 
more Wolf-Rayet stars at the expense of red-sup ergiants 
IVazauez et al.ll2007l : lEldridge. Izzard fc Toutll2008l ) 

Both can also increase the number of massive main- 
sequence stars observed. Rotation does this by increas- 
ing the amount of mixing during the main sequence and 
thus extending lifetimes and making stars of the same 
mass and luminos i ty ha ve higher surface temperatures 
l|Maeder fc Mevned |2005| ). For a detailed comparison of 
our single-star m o dels t o the Geneva rotation models see 
lEldridge fc VinkI l|2006l ). Binaries can also increase the 
number of main-sequence stars, due to secondaries ac- 
creting material during binary interactions and becoming 
more massive. However there is evide nce that in such pro- 
cesses rotation may be imp ortant (e.g. lCantiello et al.ll2007l : 
IStancliffe fc Eldridgell2009l) . The difference between the two 
processes however is that the enhancement by binaries can 
be delayed to much later times than rotation. In addition 
for rotation to have a significant effect all stars would need 
to have initial rotation velocities around 300km s~^. 

We omit rotation from our models, probing instead the 
degree to which binarity alone can explain the observed fea- 
tures of stellar populations. Discrepancies between observa- 
tion and our theoretical models may indicate the role played 
by rotation. 

2.2 Model atmospheres 

At each timestep, the properties of the stellar evolution 
model are used to select the most appropriate theoretical 
stellar atmosphere model from one of three sources. 

Firstly, for stars with hydrogen envelopes and surface 
temperatures <25kK, we use t he widely employed BaSeL 
V3.1 model atmosphere library (|Westera et al.ll2002l '). Stars 
with higher surface temperatures are treated as OB stars. 
For t hese we use the high-resolution v ersions of the mod- 
els of ISmith. Norris fc Crowtheil l|2002l fl Both libraries are 
arranged in a grid over effective temperature and surface 
gravity. We interpolate within this grid linearly to obtain 
appropriate spectra for our models. 

The most important difference between this work and 
previous analyses is o ur use of the theoretical atmospheres o f 
the Potsdam group l|Hamann. Grafener fc LiermannI |2006| ) 

^ Which can be found at 

http : //zuserver2 . star .ucl . ac.uk/ ~lj s/starburst/BM_models/ 
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for Wolf-Rayet stars (defined as having surface hydrogen 
meiss fraction X ^ 0.4 and log(rcfr/K) ^ 4). These are ad- 
vanced atmosphere models that can be related directly to 
stellar evolution models. We use not only the publicly avail- 
able models for WN stars but also a set of WC models from 
the same group that are preliminary results (W.-R. Hamann 
& A. Barniske, private communication). We have compared 
th e atmosphere models to low resol ution spectra produced 
bv lSmith, Norris fc Crowtheij (|2002l ) and find broadly simi- 
lar results to those from these more detailed models. These 
represent an important step forwards to producing a solely 
theoretical synthetic population spectrum rather than one 
based on difficult to interpret empirical observations. These 
models are on a grid of transformed radius and effective tem- 
perature which we interpolate linearly between the values of 
our stellar models. 

We use the Potsdam WR atmosphere models in the 
parameter space they cover when X ^ 0.2 and log(roff/K) 
4.45). We use the WNL models when 0.2 ^ X ^ 0.1. When 
0.1 ^ X ^ 0.01 we interpolate between the WNL models 
and the WNE models. We use the WNE models alone when 
X < 0.01. 

To determine when to switch to WC models we use the 
variable a = (xc + xo)/y where xc, xo and y are the abun- 
dance by number of carbon, oxygen and helium respectively. 
When a > 0.01 we begin to interpolate between the WNE 
and WC atmosphere models until a > 0.26. This value is 
chosen as it is the value of the composition of the atmo- 
sphere model. Above this value we use the WC atmosphere 
model alone. 

This scheme omits a subset of stars that should be in- 
cluded as WR stars, those with 4 ^ log(refr/K) ^ 4.45 
and 0.2 ^ X ^ 0.4. For these stars no model WR star 
spectra exist, so we use the corresponding BaSeL or OB 
spectra but modified to include the line luminosities of the 
Hell line at 1640A and the WR blue bump contributed by 
these stars. We use the empirical line luminosities given in 
iBrinchmann. Pettini fc CharlotI l|2008l . based on Crowther & 
Hadfield , 2006) for WNL st ars rather than the older values 
in Sch aerer fc Vaccal (|l998i V We only apply this correction 
if the star has a luminosity above log{L/ Lq) 4.9 to pre- 
vent a large contribution to these features due to lower mass 
stars at late times. By using fixed line luminosities given in 
Table [T] at our assumed luminosity limit, 10 percent of the 
total stellar emission is in these lines. This proportion would 
be larger if we included these line luminosities for less lumi- 
nous stars. In lEldridge. Izzard fc Tout! (|2008l ') it was found 
that this luminosity limit was necessary to reproduce the 
observed ratio of the number of WC to WN stars. In Table 
[2] we show the mean line luminosities predicted for all WR 
stars in our synthetic population. We find that in general we 
predict slightly lower mean line luminosities than those de- 
rived from observations, but they are similar in magnitude. 
For the WNL stars this is because the Potsdam atmosphere 
models provide smaller line luminosities than suggested by 
Table [1] The WNE line luminosities calculated from the 
Potsdam atmosphere models do agree in general with the 
line luminosities in Table [1] The largest mismatch is for the 
blue WR bump line luminosity. 

Table [5] lists the mean line luminosities at different 
metallicities. The Potsdam Wolf-Rayet atmosphere models 
only exist for solar metallicity. However it is possible to use 



Table 1. Input emission line luminosities used for WR stars with 
4 ^ log(T;,ff/K) ^ 4.45 as discussed in Section l2.2l Line strengths 
are given in units of 10^^ ergss"^. 





WNL 


WNL 


WNE 


WNE 


Metallicity 


He(II) 


Blue Bump 


He(II) 


Blue Bump 


< O.2Z0 


43 


8.3 


17 


1.7 


> O.2Z0 


247 


31 


84 


8.4 



them at other metallicities because the helium, carbon, ni- 
trogen and oxygen composition of Wolf-Rayet stars is almost 
independent of initial metallicity and is determined by the 
core nuclear burning reactions. The slight changes that do 
occur are accounted for by the lower iron opacity in the evo- 
lution models, making the surface temperature and radius 
of the modelled stars greater and smaller respectively. This 
biases the population towards earlier Wolf-Rayet spectra at 
lower metallicities. However as Table [2] shows at the lowest 
metallicities we vastly overpredict the mean line luminosities 
by using these atmo sphere models unaltered. Therefore we 
adapt the scheme of IBrinchmann. Pettini fc CharloO (|2008l ) 
in that below one-fifth solar metallicity we use the reduced 
line luminosities for WNL stars and also reduce the line 
strengths of the Hell and WR blue bump in the Potsdam 
atmosphere models by a factor of a fifth as ind icated from 
the observations of lCrowther fc Hadfieldl (|2006l ) . This leads 
to a closer match to the observed mean line luminosities. We 
were uncertain whether to use the reduced line strengths at 
a metallicity of Z = 0.004 or not. We calculated a third set 
of model spectra with the line luminosities at this metal- 
licity reduced by an intermediate value of three fifths. In 
Table [2] we see that the resulting mean line luminosities 
are only slightly less than those at the higher metallicity 
of ^ = 0.008. It is more physical to expect a gradual re- 
duction in line luminosities with metallicity rather than a 
sharp drop. Therefore in the rest of this paper we alter the 
line strengths at Z = 0.004 by three fifths in the Potsdam 
models and use the mean of the values listed in Table [1] for 
the WNL stars not covered by the Potsdam models. 

The second empirical input into the stellar spectra is 
to take ac count of Of stars. For those w e agai n use the 
method of IBrinchmann. Pettini fc Chariot (|2008l ). enhanc- 
ing the equivalent width of certain emission lines in the O 
star spectrum when the grav ity of the O star i s less than that 
of Of star s as described by Leitherer et al.l (Il999l ). As dis- 
cussed by IBrinchmann. Pettini fc CharlotI (|2008l ). Of stars 
must be accounted for to create an accurate spectrum at 
low metallicity. These are the most luminous type of O 
star and have broad emission lines, but differ from Wolf- 
Rayet stars. Therefore when the surface temperature of a 
model is greater than 33kK and the gravity is less than 
3.676 log(rofi/K) — 13.253, we supplement the Hell emis- 
sion lines at 1640A and 4686A to produce line luminosities 
of 20 and 1. 



X 10^^ ergs s ^ respectively. 



2.3 Producing a total synthetic population 
spectrum 

The procedure outlined above yields a synthetic spectrum 
appropriate to each timestep of a stellar evolution model. 
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Table 2. Mean emission line luminosities predicted by our mod- 
els for WR stars, in units of 10"^^ ergs s""*^. Metallicities marked by 
a single asterisk have had the input line luminosities reduced by a 
fifth and metallicities marked by a double asterisk have the input 
line luminosities reduced by a factor of three fifths as discussed 
in Section l2.2l The values given in parentheses are the mean line 
luminosities of WR stars where only the Potsdam model atmo- 
spheres have been used. 



WNL WNL WNE WNE 

Metallicity He(II) Blue Bump He{II) Blue Bump 



Single 

0.001* 32 (8) 3 (4) 35 0.04 

0.001 198 (103) 26 (19) 139 0.1 

0.004* 42 (20) 5 (4) 48 1 

0.004** 128 (70) 16 (10) 98 2 

0.004 214 (121) 28 (16) 148 2 

0.008 202 (105) 26 (14) 122 8 

0.020 158 (100) 22 (13) 100 17 



Binary 

0.001* 48 (14) 6 (3) 28 0.1 

0.001 218 (74) 28 (9) 102 0.4 

0.004* 46 (13) 6 (2) 26 0.6 

0.004** 148 (38) 19 (5) 58 1 

0.004 206 (63) 27 (8) 89 2 

0.008 182 (49) 24 (6) 67 4 

0.020 145 (52) 20 (8) 58 7 



To construct a synthetic spectrum for the population 
as a whole we combine the spectra for models of different 
initial mass, scaling the contribution of each by its bolo- 
metric luminosity and weighting by the timestep and initial 
mass function. Thus to scale the spectra for a specific pop- 
ulation one has simply to select the required total initial 
mass and determine the star formation history. We use two 
sets of model: single stars and binary stars. Because we are 
using detailed binary models we can follow the effects of 
binary evolution at every timestep. This is particularly ad- 
vantageous during mass-transfer events when the hydrogen 
envelope is being removed. We bin the spectra in time by 
log(age) with bins 0.1 dex wide. The youngest age we con- 
sider is 1 Myr. 

When a spectrum is added into the total synthetic pop- 
ulation we weight its contribution by a Salpeter IMF with 
a — —2.35. We calculate the total stellar mass but assuming 
minimum and maximum initial masses of 0.1 to 120 Mq. All 
the model populations presented below have an initial mass 
of 10^ Mq . For the binary populations we use the same IMF 
to determine the mass of the primary star and set the pa- 
rameters of the secondary as discussed in section 12.1.11 The 
total mass of primaries and secondaries is taken to be the 
same as for the single stars, that is 10^ Mq. 

The result of this process is a totally synthetic simu- 
lated spectrum for a young massive star population over any 
timescale and metallicity required, with any star formation 
history. 

We note that shi fting from the Salpeter I MF to the 
Kroupa IMF used by lEldridge. Izzard fc ToutI (|2008l ) has 
no measurable effect on the Hell line. Varying the IMF 
by ±0.35 can alter the Civ EW by 10 percent but only 
at ages below 2Myrs. At older ages the effect is negligible. 



A larger effect is found on the WR features in the opti- 
cal spectrum. We find that the optical WR features here 
can vary by up to 5 percent as we alter the number of 
low mass main sequence stars contributing to the optical 
continuum relative to the number of WR stars. Such small 
ch anges in the IMF also have littl e effect on the ratios given 
in lEldridge. Izzard fc ToutI (|2008l ). 



2.4 Taking account of nebular emission and other 
details 

One final detail in our spectral synthesis is to include 
the contribution from nebular emission. In star-forming 
galaxies, interstellar gas is ionised by the stellar contin- 
uum emitted blueward of 912A, and upon recombina- 
tion it emits a nebular continuum. Neglecting this emis- 
sion would lead to an incorrect estimate of the equivalent 
widths of emission lines and incorrect broad-band co lours 
(jZackrisson. Bergvall fc Leitetllioosi : iMoUa et al.ll2009l '). We 
use the program Cloudy to produce a detailed model of 
the output spectrum from our stellar spectra. We give the 
details of the Cloudy models we use below. The model out- 
put is sensitive to the chosen inner radius and composition 
of the gas used in the code. The details of our illustrative 
nebular emission model are as follows: 

• metals (Z/0.02) linear 

• element scale factor hydrogen ((0.75 — 2.5.Z)/0.7) 
linear 

• element scale factor helium ((0.25 + 1.5Z)/0.28) 
linear 

• hden 2 log constant density 

• covering factor 1.0 linear 

• filling factor 1.0 linear 

• sphere 

• radius 1 . log parsec 

• iterate 

• set temperature floor 1000 

• stop temperature lOOK 

• stop efrac -2 

The details can be altered to reproduce a more accurate 
model spectrum however we only use this simple scheme here 
to demonstrate how the synthetic spectrum is affected by 
nebular emission in Section 3.1. We do not attempt to model 
the nebula emission lines when we compare our models to 
observed spectra in Section 3.3 below. To match the neb- 
ula emission lines requires altering the input of our Cloudy 
model and only tell us the composition of the nebula gas 
while in this paper we are primarily interested in the stel- 
lar emission. We also omit dust from our Cloudy models; 
we note that we are only concerned with young stellar pop- 
ulations and with the ratio between line emission and the 
continuum. Assuming a uniform dust geometry, these will 
be suppressed equally leaving the line equivalent width s un- 
affected. However as discussed by I Chariot fc Fall (|200(]| ) and 
IConrov. Gunn fc" White! I 20081 ) for galaxies made up of mul- 
tiple stellar population the situation may be more compli- 
cated as the birth clouds may disperse on average after 10 
million years and so older populations may be attenuated 
less by dust than younger populations and therefore domi- 
nate the observed spectrum. 



6 Eldridge & Stanway 



The theoretical model spectra in their raw state repre- 
sent a population in which all the stars are static, lying at 
precisely the same redshift or recession velocity, and makes 
no account for the velocity dispersion of stars within the ob- 
served galaxies. To account for broadening due to the veloc- 
ity dispersion of stars within a typical star-bursting galaxy 
(< lOOkms"^), we convolve our final spectra with a boxcar 
function of width lA. For strongly star-forming galaxies at 
higher redshifts a higher velocity dispersion might be appro- 
priate. 



3 VERIFICATION AND LOCAL 
APPLICATION 

Any new model must be tested against observations. Here we 
compare our models to three different sets of observations on 
sites of star- formation expected to have WR stars present. 



3.1 Unresolved young massive star clusters in 
nearby galaxies 

Before we compare our models to distant galaxies it is sen- 
sible to compare them to nearby unresolved stellar popu- 
lations. Therefore we compare our models to a large set 
of observations of young massive star clusters compiled by 
iLarsen fc Richtlerl (|l999l ). They observed a number of star 
clusters in nearby spiral galaxies with broad-band photome- 
try. Since the host galaxies were spiral we assume the metal- 
licity of the clusters is broadly solar and compare the ob- 
servations to colours from our models for a cluster with a 
mass of 10^ Mq. For comparison we have created a similar 
track from starburst99 with the same total stellar mass and 
initial mass distribution. T his track was calculated using 
'Starburst99 for Windows' (|Leitherer fc ChenI [20091 ') using 
the same IMF and total mass as for our models, the stan- 
dard Geneva solar metallicity stellar evolution tracks, the 
Lejeune stellar atmosphere models and the remaining op- 
tions at their default values. We plot the results in Figure 
[1] Our tracks are shorter than the starburst99 results as we 
terminate our simulations at 40 Myrs since we do not cur- 
rently include stars with initial masses below 5 Mq which 
become important at times later than this. 

From Figure [T] we see that the evolution of colours is 
broadly consistent between models and observations. How- 
ever there are a number of important differences between 
the model tracks. These differences are primarily due to 
the different stellar models employed by starburst99 and 
the models presented here. Our model tracks tend to pass 
through regions of the diagram that contain more observed 
clusters although the observations have a large scatter, simi- 
lar to the random error of the observations. In B-V the star- 
burst99 models tend to have redder colours than our models 
at around 10 Myrs. However we see the greatest differences 
between the different model tracks in V-I. Our single star 
population tracks extend much further into the red than our 
binary models because binary models reduce the number of 
red supergiants and therefore reduce their contribution in 
the I band. However, the overproduction of red supergiants 
is a feature of all models containing only single stars (with 
the main difference being how those supergiants are charac- 



terised in the model population) and emphasises the need 
for binary population models. 

Our binary models also tend to be bluer in B-V than 
our single star models, this is due to the increase in the 
number of main-sequence stars at late times due to secon- 
daries accreting mass in binary interactions. We note that 
all the models plotted, those presented in this paper and 
similar models from starburst99, deviate from the locus of 
observed clusters in the B-V vs U-B colour plane at early 
times (< 5Myr), most likely due to variations in the metal- 
licity of the clusters away from the Solar composition of our 
tracks. This deviation may also be due to relatively simple 
approximations for the strong nebular continuum contribu- 
tion at these times. Models omitting this nebular contribu- 
tion, while unphysical, can provide a better fit to the data in 
this region, suggesting that more detailed modelling of the 
nebular emission may be necessary. 

There are a number of factors that we have not included 
in our model tracks presented here. For example we have 
not attempted to fit absorption from dust in the HII region 
Cloudy model. This is because the starburst99 model track 
does not include dust, only nebula continuum emission. In 
each of the panels we indicate the reddening direction for 
an Av of 0.5. We see that in some of the scatter of the 
observations could be explained by line of sight dust. 

Given these caveats, we nonetheless find that binary 
populations have very different colours to a single star pop- 
ulation at certain ages, and that the observational data at 
these points is often more consistent with the binary mod- 
els with few sources showing V-I colours, for example, as 
extreme as those predicted for single star populations. To 
achieve this difference between the single and binary mod- 
els, binaries with orbital separations between 100 to 1000 Kq 
must be included. Wider binaries do not interact so produce 
results little different from those from single stars. Tighter 
binaries tend to experience mergers and so, while evolving 
as binaries for some of their lives, they eventually become 
single stars. 

3.2 Unresolved stellar populations from SDSS 

iBrinchmann. Kunth fc Purred (|2008l ) recently presented a 
study of a selection of SDSS galaxies showing evidence for 
ongoing massive star formation. They searched the SDSS 
DR6 archival spectra to identify those with Wolf-Rayet fea- 
tures in the optical. The two features used for identification 
are known as the blue and red Wolf-Rayet bumps at approx- 
imately 4700A and 5800A. These star-forming galaxies host 
massive stellar populations similar to those incorporated in 
our models, and hence provide a good experimental verifi- 
cation of the predictive power of our models. In order to 
construct an appropriate stellar population in our synthetic 
spectrum, we assume an instantaneous burst of star forma- 
tion and consider its evolution with tim^fl- We calculate the 
strength of Wolf-Rayet features both without and with a 
nebula emission model created with Cloudy as detailed in 
Section 12.41 Inclusion of nebular continuum emission leads 

^ Note that this d iffers from our earlier work 
llEIdridge. Izzard fc Tou^ [ 20081) . which considered continu- 
ous star formation. 
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Figure 1. The theoretical broad-band colours for a stellar popula t ion w ith a total mass of IO^A/q compared to observations of unresolved 
young massive clusters in local galaxies from iLarsen fc Richtled l|l999l V The theoretical models include a prediction from starburst99, 
our single-star population and our binary population all including nebular emission at solar metallicity. The change in colours produced 
by extinction of = 0.5 is indicated by a vector in each case. Tickmarks on the colour tracks indicate age increments of 0.1 dex, from 
IMyr at the bluest B-V colours to 40 Myr. 



to only a slight decrease in the maximum EW allowed by up 
to 10 percent. 

Here we calculate predicted values for the strength 
(equivalent width; EW) of Wolf Rayet features as defined 
by Brinchmann et al, removing the contribution of nebular 
lines to allow fair comparison with the data. As illustrated 
by Figures [2] and [3] in general our models reproduce the 
range of observed equivalent widths of both features with 
the exception of a few extreme cases. 

For a simple check of our nebula model, Figure|3]shows a 
comparison between the equivalent widths of the Blue WR 
bump to a nearby nebula emission line, H-/3. Comparing 
the ratio of these EWs will indicate whether we are esti- 
mating the strength of the nebula emission lines correctly 
relative to the stellar spectrum. We see that our model pre- 
dictions of this ratio agree with the range of ratios from ob- 
serve d galaxies and previous estimates from the binary mod- 
els of IVan Bever fc Vanbevereiil l|2003l ). Binary populations 
at metallicities below solar metallicity produce the highest 
ratio we discuss why this is below. We note our predicted 
ratios are a lower limit as the ratio we predict can be varied 
by changing the details of our Cloudy model, especially the 
covering factor. Decreasing this lets more ionising photons 
to escape and therefore decreases the nebula emission fea- 



tures such as the H-/3 line luminosity without significantly 
affecting the blue WR bump EW which is determined by the 
stellar population. Thus the vertical spread has degeneracy 
between age and covering factor. 

In Figure Owe show the strength of the key Wolf-Rayet 
emission features in our model spectra as a function of the 
age of the stellar population and its metallicity. We see that 
the spread of values in Figures and U is due to a range of 
ages in the stellar population. Emission features peak at ages 
of a few million years for single stars and can be extended 
to much later times, up to 10 million years, by the inclusion 
of binaries. We note that similar diagrams with a nebula 
continuum component would reduce these observed EWs by 
a small amount. In Figure |S] we show how the strength of 
H-/3 varies with age. We see that binary models exhibit H-/3 
to later ages. This is why at solar metallicity in Figure |3] at 
Solar metallicity single stars have a greater maximum ratio 
than binary stars, a smaller H-/3 EW gives a larger ratio 
rather than great blue WR bump EW. 

In Figure[5]there are large peaks in the strength of Hell 
and the blue WR bump at around 10 million years. These 
peaks are due to WNL stars formed from binary evolution. If 
these stars were in binaries they would be red supergiants. 
Their WNL phase has a long lifetime due to weak stellar 
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Figure 2. The EW of the blue WR bump measured 
in SPSS galaxies and corrected for nebular emission by 
iBrinchmann. Kunth &: Durretl ||2008|) compared to predicted EW 
from our models without nebula emission. Small black symbols 
indicate observations. The models are represented by + for single- 
star populations and X for binary populations. The colour coding 
for the models are red: twice solar, yellow: solar, green: two-fifths 
solar, cyan: one-fifth solar and blue: one-twentieth solar. The ver- 
tical spread is due to the change of the EW with age of the stellar 
population. 



winds at lower metallicities, so it takes longer for stellar 
winds and binary evolution to remove their hydrogen en- 
velopes than at higher metallicity. Therefore their contribu- 
tion to the composite spectrum is greater. Also as seen in 
Figure |5] the H-/3 EW is small therefore artificially boosting 
the ratio of the two EW in Figure |4] If we only consider 
the ratio at ages less than lO'''* years we do not achieve the 
highest ratios in Figure |4] and the ratio decreases with the 
same trend indicated by observations. 

The reasons why such low metallicity systems with 
large ratios ar e not observed are related to the work of 
ICharlot fc Falll (|2000l ) suggesting the gas in a massive clus- 
ter disperses by around 10 million years. With little gas sur- 
rounding the stars most ionising photons will escape produc- 
ing little observation H-/3 emission, effectively decreasing the 
covering factor. In the sample of SDSS galaxies we use only 
a few galaxies have H-/3 EW less than lOA. If the cover- 
ing factor decr e ases w ith ages as suggested by the work of 
ICharlot fc Falil (|2000l ). our predicted H-/3 EW will drop be- 
low this value. Therefore it becomes less likely to observe the 
large ratios of blue WR bump EW to H-/3 EW we predict 
in Figure |4] 



3.3 Comparison to the spectra of Tol A &: B 

A further comparison between our models and the inte- 
grated spectra of young stellar populations can be made 
using star formi ng regions with good spectral a nd photo- 
metric coverage. ISidoli. Smith fc Crowtheij (|2006l ) observed 
the massive stellar population in the giant HII region Tol89 
in NGC5398. The published spectra of these regions span 
from the UV (useful for comparison to high redshift sources) 
to the optical (easily observed in more local galaxies). 
ISidoli. Smith fc Crowtherl (|2006l ) attempt to fit the observed 
Wolf-Rayet lines with a starburst99 model. They find that 
for the two knots of star formation, A & B, the ages are 



Figure 3. Similar to Figure [2] but for the EW of the red Wolf- 
Rayet bump. 
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Figure 4. Similar to Figure [2l and [3l but now for the ratio of the 
Blue WR bump to the H-/3 emission line. The H-/3 EW is esti- 
mated from our Cloudy and is an upper limit as the H-/3 EW can 
be reduced by decreasing the covering factor in the model. There- 
fore the ratios here are a lower limit as decreasing the covering 
factor will increase the ratio. 

4. 5 ±1.0 and < 3 Myrs respectively with an LMC-like metal- 
licity, roughly two-fifths solar, preferred. 

We attempt to fit the same spectra with our models. 
To do this we first measure the equivalent widths of the ul- 
traviolet Hell A1640A emission line and the Red and Blue 
Wolf-Rayet bumps in an identical manner on both the ob- 
served and model spectra. The Hell EW is calculated by nor- 
m alising the spectr a using the continuum points identified 
bv lRix et ai] l|2004l 'l and estimating the EW over the wave- 
length range 1635 to 1650A. We calculate EWs for the WR 
bumps as in Section 13.21 but studying the spectrum to re- 
move any (relatively narrow) nebula lines. The UV spectrum 
of Tol A is a stack of the spectra from 4 individual knots of 
star formation. We find one of these knots, A4 as identified 
bylSidoh, Smith & Crowthcr (2006) has a very narrow neb- 
ula emission of Hell at 1640A. This dominates the EW of 
the A4 knot and has an equivalent width of 0.2A in the Tol 
A stack. Nebula Hell emission indicates a very hard ionis- 
ing radiation field in this cluster. It has been suggested such 
lines indicate the presence of population HI stars ( Scha crerl 
l2003h . However in this case the hard ionising field is more 
likely to come from some other source given the relatively 
high metallicity and zero redshift. For example that knot is 
old enough that there could be an accreting black hole or 
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Figure 5. The EW of the main Wolf-Rayet spectral features 
measured from our model spectra as a function of age and metal- 
licity. Solid lines are for single stars models, the dashed lines are 
for binary star models. The colour scheme for the metallicities is 
given in Figure [2] The horizontal black lines represent the mea- 
surement of these lines for Tol A & B. The solid line is for Tol 
A with the uncertainty shown by the striped region. The dashed 
line is for Tol B with the grey region showing the uncertainty. 
These synthetic spectra do not include a nebular emission model. 



neutron star present producing harder radiation. We have 
not yet attempted to include such emission in our models 
yet. Therefore we have calculated the EW of the Hell line 
of Tol A ignoring the A4 knot. 

We plot our model line strengths against age in Figure 
[5] compared to the observed EW of features in the Tol A & 
B regions. We find that single star and binary models are 
both able to reproduce the observed EW in these features 
simultaneously at around 10^'^ years. We also show how 
show how the H-/3 EW varies with age in Figure [6] While 
we find that the EW broadly agrees with those observed for 



2.5 




bg(Age/yrs) 

Figure 6. Similar to Figure [5] but now for the EW of the H-/3 
emission line. The H-/3 EW is estimated from our Cloudy models 
and is an upper limit as the H-/3 EW can be reduced by decreasing 
the covering factor in the model. While we do not use this EW to 
determine an age for Tol-A & B we note that the EW from the 
observed spectra are log(EW(H— /3)) = 1.6 and 2.2 respectively 
which broadly agree with our predicted H-/3 EW at the ages we 
derive from the stellar emission lines. 



Tol A & B we do not use them to derive the ages of the 
stellar population. 

The best fit metallicity is an LMC-like metallicity of 
O.4Z0. The age of Tol A is constrained to be from 10®'^ and 
lO'' '^ years. At ages outside this range the stellar features 
are too weak to agree with the observed values, or for the 
binary models, too strong. For Tol B the age is constrained 
to be less than 10®'^ years. 

To verify the best fit ages and metallicities we show 
a more detailed comparison between these models and the 
observed spectra for Tol A and B in Figures [7] and |8] re- 
spectively, considering the three key WR diagnostic features 
and also the ultraviolet CIV absorption/emission feature at 
1550A, which is sensitive to the presence of massive stars 
and in particular winds driven from O-stars. For Tol A, we 
find reasonable visual fits between 10^'^ and 10^'® years at 
both SMC (O.2Z0) and LMC metallicities. The CIV line re- 
quires younger ages for the system than the WR features 
provide. For Tol B we find any model with negligible Hell 
emission provides a reasonable fit to the data. However mod- 
els with ages of 10®'^ years and older tend to overpredict the 
CIV emission line relative the observed spectrum. Also for 
the blue and red bumps there are no distinctive broad emis- 
sion lines and the spectrum is dominated by narrow emission 
lines. This again suggests a young population. 

Our derived ages, both from matching the spectra by 
eye and from the line equivalent widths, are 4 Myrs for Tol 
A and < 2.5Myrs for Tol B. Th ese agree with the ages of 
ISidoU, Smith fc Crowtheil (|2006l ). However our ages are on 
the young side of the ages they derived. 

By eye the model agreement is good, especially in light 
of the general models employed, which are not fine-tuned to 
any great degree. The line profiles and strengths of ultra- 
violet spectral features are generally less well fit than the 
optical features. There are several ways in which our mod- 
els might be fine tuned to fit the observational data in this 
case. Firstly, it may be possible to refine our model for neb- 
ular continuum contribution to the integrated spectrum. At 
present, we consider a basic model as discussed above and 
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parameters could be adjusted to gain a better fit for this 
star-forming region. Secondly, it may be possible to include 
the effects of material external to the system in addition to 
gas and dust local to the star-forming region which is mod- 
elled using Cloudy. Thirdly, we only use a sparse grid of 
O-star models rather than a dense grid of atmos pheres and 
it ma y be possible to improve upon this in future l|Rix et al.l 
I2OO4I ). Fourthly, we have assumed an instantaneous burst. 
A better fit might be achieved by allowing multiple star for- 
mation episodes with slightly different ages. Fifthly we have 
not considered the effects of stellar rotation on the evolu- 
tion or the spectra. This will have two effects on our results. 
It would extend the main-sequence life-time and would also 
rotationally broaden the CIV line further, particularly if the 
stars rotate at velocities above 200km s^^. This would po- 
tentially improve the fit between models and observed data. 
We speculate that study of the CIV profile could potentially 
give a way to evaluate the importance of rotation versus bi- 
narity in stellar populations. Finally, some WNL stars have 
CIV in emission. We have not included any such CIV emis- 
sion in the WNL stars that are not represented by the Pots- 
dam model atmospheres. Including an empirical correction 
for such emission in these models may broaden the emission 
component of the CIV line and simultaneously decrease the 
absorption component. However we do not choose to intro- 
duce arbitrary line emission here. 



4 DISCUSSION AND CONCLUSIONS 

In this paper we set out to describe the construction of model 
stellar populations incorporating massive stars and massive 
stellar binaries, and then the synthesis of spectra for this 
population. We have compared our model spectra to ob- 
servations of comparable unresolved stellar populations and 
found agreement is fair over a large range of wavelengths, 
from the UV to the near infra-red. 

In general our population synthesis including binary 
models predict that such systems have less emission at long 
wavelengths around the I band. This is because binary in- 
teractions remove the hydrogen envelope of some red super- 
giants to form Wolf-Rayet stars. These WR stars then lead 
to more blue colours in B-V and V-I broad-band colours 
and a larger UV flux. The latter increases the timespan over 
which nebula emission is important to the evolution of stellar 
populations from 6 Myrs to 20 Myrs. Another binary effect 
is that it is more likely that strong WR emission lines are 
observed since binary interactions tend to spread WR emis- 
sion features over a longer timespan. For single stars the WR 
stars all exist over a short timespan so WR features are only 
present for a short period of time. This suggests that ages 
derived from Wolf-Rayet features to date may have been 
systematically underestimated. 

It may be somewhat surprising that by including bi- 
nary evolution we find no great difference from predic- 
tions from starburst99, a code based on single-star evolu- 
tion alone. However the stellar evolution models used in 
the starburst 99 model present in Section TS. II are nearly two 
decades old (|Schaller et al.lll992l ). The important difference 
between these models and those used here is that the mass- 
loss rates have substantially decreased for OB-stars and WR 
stars. Therefore with incorrect single-star mass-loss rates 



the older models reproduce the observed stellar population. 
Our model binary population therefore should broadly agree 
with this older code, but now we model the same magnitude 
of mass-loss as a combination of lower single-star mass-loss 
rates enhanced for some stars by binary interactions which 
is a physically reasonable model. 

We demonstrate that our code produces a good fit to 
the observational data of local stellar populations in which 
massive stars are important. However, further refinements 
are possible and additional verification data on local star- 
forming regions would be welcome. Nonetheless, the stellar 
models and synthesis code presented here may now be used 
as a tool to study stellar populations in a range of different 
observational domains and to derive their physical parame- 
ters, as demonstrated by their application to extragalactic 
star- forming regions Tol A and B. 
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Figure 8. Similar to Figure [7] our model spectra for Tol B, the left panel is for a model with SMC-like metallicity and an age of 
10^'^ years and the right panel is for a model with LMC-like metallicity at an age of 10^' years. 
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